clear;clc;
rp = 100e-9;
x = [4,6,8,10,12,13,14,16,20];
l = [402.1087, 429.5096, 458.2806, 485.6815, 514.4525, 528.1529, 541.8534, 569.2543, 621.3161] * 1e-9;
r = [40.4, 40.1, 39.7, 39.4, 39.2, 39.1, 38.9, 38.7, 38.4] * 1e-9;
k = [5.0519E-18, 8.3021E-18, 1.0335E-17, 1.1778E-17, 1.2549E-17, 1.2907E-17, 1.3031E-17, 1.3244E-17, 1.3523E-17];
c = (l ./ (l-2*rp) + l ./ (sqrt(2)*l-2*rp) * 4 + l ./ (sqrt(3)*l-2*rp) * 4) / 9;
k0 = (r .^ 4) ./ (l .^ 2) .* c;
y = k ./ k0;
cla;
plot(x, y, 'kx', 'MarkerSize', 10, 'Linewidth', 2);
hold on;
p = polyfit(x, y, 1);
x3 = 0:26;
y3 = polyval(p, x3);
corr = corrcoef(x3, y3);
plot(x3, y3, 'r', 'Linewidth', 2);
axis([0, 20, 0, 3]);
title('Coordination Number Function');
xlabel('Coordination Number \zeta');
ylabel('f_0(\zeta)');
text(11, 1, ['f_0(\zeta) = ', num2str(p(1), '%.2f'), '(\zeta - ', num2str(-p(2)/p(1), '%.2f'), ')'], 'Fontsize', 14);
text(11, 0.7, ['R^2 = ', num2str(corr(1, 2), '%.4f')], 'Fontsize', 14);
legend('Simulation', 'Linear Continuation');